rm(list = ls())
library(dplyr)
library(wesanderson)
library(GillespieSSA)
library(tidyverse)
Information regarding births, deaths and population size were obtain from a study conducted by Headland et al., (2011). Authors conducted a census-like survey of the
agta_demo <- read.csv("AgtaPopDynamics_Headland2007.csv")
ggplot(agta_demo, aes(x=Year)) +
geom_line(aes(y=PopSize), colour = wes_palettes$Darjeeling1[1]) +
geom_line(aes(y=Births), colour = wes_palettes$Darjeeling1[2]) +
geom_line(aes(y=Deaths), colour = wes_palettes$Darjeeling1[3]) +
theme_bw()
hist(agta_demo$PopSize)
summary(agta_demo$PopSize)
Min. 1st Qu. Median Mean 3rd Qu. Max.
133.0 177.0 213.0 211.4 228.0 295.0
hist(agta_demo$Births)
summary(agta_demo$Births)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
4.00 8.00 10.00 10.33 12.25 15.00 1
hist(agta_demo$Deaths)
summary(agta_demo$Deaths)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2.000 5.000 7.000 7.683 10.000 23.000 1
agta_demo <- agta_demo %>%
mutate(Birth_rate = Births/(Female*2),
Birth_rate_daily = (1 + Birth_rate) ^ (1/365) - 1,
Death_rate = (Deaths/PopSize),
Death_rate_daily = (1 + Death_rate) ^ (1/365) - 1,
PopChange = (diff = PopSize - lag(PopSize, default = first(PopSize))),
PopChange_rate = abs(PopChange)/PopSize,
PopChange_rate_daily = (1 + PopChange_rate) ^ (1/365) - 1)
head(agta_demo)
hist(agta_demo$Birth_rate)
hist(agta_demo$Death_rate)
ggplot(agta_demo, aes(x=Year)) +
geom_line(aes(y=Birth_rate), colour = wes_palettes$Darjeeling1[2]) +
geom_line(aes(y=Death_rate), colour = wes_palettes$Darjeeling1[3]) +
theme_bw()
demo_sum <- agta_demo %>%
select(PopSize, Birth_rate, Birth_rate_daily, Death_rate, Death_rate_daily, PopChange_rate, PopChange_rate_daily) %>%
summarise(across(
.cols = is.numeric,
.fns = list(Mean = mean, SD = sd), na.rm = TRUE,
.names = "{col}_{fn}"
))
demo_sum
demo_sum <- as.list(demo_sum)
# Import Camp data from Mark Dyble
camps.data <- read_csv("camps.csv")
New names:Rows: 15 Columns: 9── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): camp_name
dbl (8): ...1, camp_total, camp_adult_men, camp_adult_women, camp_all_r, camp_adult_r, forage, turnover
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(camps.data)
# Explore camp size
hist(camps.data$camp_total)
camp.size <- camps.data %>%
summarise(mean = mean(camp_total),
sd = sd(camp_total),
min = min(camp_total),
max = max(camp_total),
var = var(camp_total))
camp.size
# Define Paramenters
N <- sample(camps.data$camp_total, 1) # Population size
initial_infected <- 1 # Initial infected
simName <- "SEIRS model" # Simulation name
tf <- 365*3
#Collect parameters
parms <- list(
beta = 0.6,
sigma = 0.175, # E to I rate
gamma = 0.2, # I to R rate
omega = 1/180, # R to S rate
mu = demo_sum$Birth_rate_daily_Mean, # Birth/death rate per person per day
alpha = 1/1000)
#Create the named initial state vector for the U-patch system.
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
# Define the state change matrix for a single patch
nu <- matrix(c( -1, 0, 0, +1, +1, -1, 0, 0, 0, 0, # S
+1, -1, 0, 0, 0, 0, -1, 0, 0, 0, # E
0, +1, -1, 0, 0, 0, 0, -1, 0, -1, # I
0, 0, +1, -1, 0, 0, 0, 0, -1, 0, # R
0, 0, 0, 0, +1, -1 ,-1, -1, -1, -1), # N
nrow=5,byrow=TRUE)
# Define propensity functions
a <-c(
paste0("(beta*I/N)*S"), # Infection
paste0("sigma*E"), # Becomes infecious
paste0("gamma*I"), # Recovery from infection
paste0("omega*R"), # Loss of immunity
paste0("mu*N"), # Births
paste0("mu*S"), # Deaths (S)
paste0("mu*E"), # Deaths (E)
paste0("mu*I"), # Deaths (I)
paste0("mu*R"), # Deaths (R)
paste0("alpha*I") # Deaths from infection
)
Define functions to calculate R0 and expected number of susceptibles at equilibrium, and critical community size (Diekmann et al., 2012).
epsilon^2
[1] 4.610078e-07
CCS
[1] 4882736
plot_data %>%
filter(state == "I") %>%
slice_max(count)
ggsave(filename = "single_plot.pdf",
plot = single_plot,
device = "pdf",
width = 7,
height = 5,
path = "/Users/matthewhoyle/Github_R_projects/Hunter_Gatherer_models")
```{r}
Error: attempt to use zero-length variable name
# Summary table of endpoint data
sim_output <- sim_output %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output
# Make Summary Table of output
sim_summary <- sim_output %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100,
mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/180)
sim_summary
#Collect parameters
parms_0 <- parms
parms_0$omega <- 0
# Run simulations with the Direct method
set.seed(4)
out_0 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_0,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_0 <- out_0$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_0 <- ggplot(data = plot_data_0, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_0
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_0 <- list()
sim_list_0 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_0 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_0,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_0 <- out_100_0$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_0[[i]] <- sim_data_0
}
sim_output_0 <- bind_rows(sim_list_0)
# Summary table of endpoint data
sim_output_0 <- sim_output_0 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_0
# Make Summary Table of output
sim_summary_0 <- sim_output_0 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100) %>%
mutate(omega = 0)
sim_summary_0
#Collect parameters
parms_1 <- parms
parms_1$omega <- 1
# Run simulations with the Direct method
set.seed(4)
out_1 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_1,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_1 <- out_1$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_1 <- ggplot(data = plot_data_1, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_1
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_1 <- list()
sim_list_1 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_1 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_1,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_1 <- out_100_1$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_1[[i]] <- sim_data_1
}
sim_output_1 <- bind_rows(sim_list_1)
# Summary table of endpoint data
sim_output_1 <- sim_output_1 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_1
# Make Summary Table of output
sim_summary_1 <- sim_output_1 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1)
sim_summary_1
#Collect parameters
parms_3 <- parms
parms_3$omega <- 1/3
# Run simulations with the Direct method
set.seed(4)
out_3 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_3,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_3 <- out_3$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_3 <- ggplot(data = plot_data_3, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_3
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_3 <- list()
sim_list_3 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_3 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_3,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_3 <- out_100_3$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_3[[i]] <- sim_data_3
}
sim_output_3 <- bind_rows(sim_list_3)
# Summary table of endpoint data
sim_output_3 <- sim_output_3 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_3
# Make Summary Table of output
sim_summary_3 <- sim_output_3 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/3)
sim_summary_3
#Collect parameters
parms_7 <- parms
parms_7$omega <- 1/7
# Run simulations with the Direct method
set.seed(4)
out_7 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_7,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_7 <- out_7$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_7 <- ggplot(data = plot_data_7, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_7
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_7 <- list()
sim_list_7 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_7 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_7,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_7 <- out_100_7$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_7[[i]] <- sim_data_7
}
sim_output_7 <- bind_rows(sim_list_7)
# Summary table of endpoint data
sim_output_7 <- sim_output_7 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_7
# Make Summary Table of output
sim_summary_7 <- sim_output_7 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/7)
sim_summary_7
#Collect parameters
parms_10 <- parms
parms_10$omega <- 1/10
# Run simulations with the Direct method
set.seed(4)
out_10 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_10,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_10 <- out_10$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_10 <- ggplot(data = plot_data_10, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_10
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_10 <- list()
sim_list_10 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_10 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_10,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_10 <- out_100_10$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_10[[i]] <- sim_data_10
}
sim_output_10 <- bind_rows(sim_list_10)
# Summary table of endpoint data
sim_output_10 <- sim_output_10 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_10
# Make Summary Table of output
sim_summary_10 <- sim_output_10 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/10)
sim_summary_10
#Collect parameters
parms_20 <- parms
parms_20$omega <- 1/20
# Run simulations with the Direct method
set.seed(4)
out_20 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_20,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_20 <- out_20$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_20 <- ggplot(data = plot_data_20, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_20
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_20 <- list()
sim_list_20 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_20 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_20,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_20 <- out_100_20$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_20[[i]] <- sim_data_20
}
sim_output_20 <- bind_rows(sim_list_20)
# Summary table of endpoint data
sim_output_20 <- sim_output_20 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_20
# Make Summary Table of output
sim_summary_20 <- sim_output_20 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/20)
sim_summary_20
#Collect parameters
parms_30 <- parms
parms_30$omega <- 1/30
# Run simulations with the Direct method
set.seed(4)
out_30 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_30,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_30 <- out_30$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_30 <- ggplot(data = plot_data_30, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_30
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_30 <- list()
sim_list_30 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_30 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_30,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_30 <- out_100_30$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_30[[i]] <- sim_data_30
}
sim_output_30 <- bind_rows(sim_list_30)
# Summary table of endpoint data
sim_output_30 <- sim_output_30 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_30
# Make Summary Table of output
sim_summary_30 <- sim_output_30 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/30)
sim_summary_30
#Collect parameters
parms_40 <- parms
parms_40$omega <- 1/40
# Run simulations with the Direct method
set.seed(4)
out_40 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_40,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_40 <- out_40$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_40 <- ggplot(data = plot_data_40, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_40
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_40 <- list()
sim_list_40 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_40 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_40,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_40 <- out_100_40$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_40[[i]] <- sim_data_40
}
sim_output_40 <- bind_rows(sim_list_40)
# Summary table of endpoint data
sim_output_40 <- sim_output_40 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_40
# Make Summary Table of output
sim_summary_40 <- sim_output_40 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/40)
sim_summary_40
#Collect parameters
parms_50 <- parms
parms_50$omega <- 1/50
# Run simulations with the Direct method
set.seed(4)
out_50 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_50,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_50 <- out_50$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_50 <- ggplot(data = plot_data_50, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_50
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_50 <- list()
sim_list_50 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_50 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_50,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_50 <- out_100_50$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_50[[i]] <- sim_data_50
}
sim_output_50 <- bind_rows(sim_list_50)
# Summary table of endpoint data
sim_output_50 <- sim_output_50 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_50
# Make Summary Table of output
sim_summary_50 <- sim_output_50 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/50)
sim_summary_50
#Collect parameters
parms_60 <- parms
parms_60$omega <- 1/60
# Run simulations with the Direct method
set.seed(4)
out_60 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_60,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_60 <- out_60$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_60 <- ggplot(data = plot_data_60, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_60
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_60 <- list()
sim_list_60 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_60 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_60,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_60 <- out_100_60$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_60[[i]] <- sim_data_60
}
sim_output_60 <- bind_rows(sim_list_60)
# Summary table of endpoint data
sim_output_60 <- sim_output_60 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_60
# Make Summary Table of output
sim_summary_60 <- sim_output_60 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/60)
sim_summary_60
#Collect parameters
parms_70 <- parms
parms_70$omega <- 1/70
# Run simulations with the Direct method
set.seed(4)
out_70 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_70,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_70 <- out_70$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_70 <- ggplot(data = plot_data_70, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_70
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_70 <- list()
sim_list_70 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_70 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_70,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_70 <- out_100_70$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_70[[i]] <- sim_data_70
}
sim_output_70 <- bind_rows(sim_list_70)
# Summary table of endpoint data
sim_output_70 <- sim_output_70 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_70
# Make Summary Table of output
sim_summary_70 <- sim_output_70 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist)/num_sims)*100,
mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/70)
sim_summary_70
#Collect parameters
parms_80 <- parms
parms_80$omega <- 1/80
# Run simulations with the Direct method
set.seed(4)
out_80 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_80,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_80 <- out_80$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_80 <- ggplot(data = plot_data_80, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_80
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_80 <- list()
sim_list_80 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_80 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_80,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_80 <- out_100_80$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_80[[i]] <- sim_data_80
}
sim_output_80 <- bind_rows(sim_list_80)
# Summary table of endpoint data
sim_output_80 <- sim_output_80 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_80
# Make Summary Table of output
sim_summary_80 <- sim_output_80 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/80)
sim_summary_80
#Collect parameters
parms_150 <- parms
parms_150$omega <- 1/150
# Run simulations with the Direct method
set.seed(4)
out_150 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_150,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_150 <- out_150$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_150 <- ggplot(data = plot_data_150, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_150
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_150 <- list()
sim_list_150 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_150 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_150,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_150 <- out_100_150$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_150[[i]] <- sim_data_150
}
sim_output_150 <- bind_rows(sim_list_150)
# Summary table of endpoint data
sim_output_150 <- sim_output_150 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_150
# Make Summary Table of output
sim_summary_150 <- sim_output_150 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/150)
sim_summary_150
#Collect parameters
parms_365 <- parms
parms_365$omega <- 1/365
# Run simulations with the Direct method
set.seed(4)
out_365 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_365,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_365 <- out_365$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_365 <- ggplot(data = plot_data_365, aes(x=t, y=count, colour=state))+
geom_line()+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_365
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_365 <- list()
sim_list_365 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
N <- sample(camps.data$camp_total, 1) # Sample different patch sizes for each sim
x0 <- c(N - initial_infected, initial_infected, 0, 0, N)
names(x0) <- c("S","E","I", "R", "N")
out_100_365 <- ssa(
x0 = x0,
a = a,
nu = nu,
parms = parms_365,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_365 <- out_100_365$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "state", values_to = "count") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, I, N, persist)
sim_list_365[[i]] <- sim_data_365
}
sim_output_365 <- bind_rows(sim_list_365)
# Summary table of endpoint data
sim_output_365 <- sim_output_365 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_365
# Make Summary Table of output
sim_summary_365 <- sim_output_365 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/365)
sim_summary_365
waning_results_single <- sim_summary %>%
bind_rows(sim_summary_1) %>%
bind_rows(sim_summary_3) %>%
bind_rows(sim_summary_7) %>%
bind_rows(sim_summary_10) %>%
bind_rows(sim_summary_20) %>%
bind_rows(sim_summary_30) %>%
bind_rows(sim_summary_40) %>%
bind_rows(sim_summary_50) %>%
bind_rows(sim_summary_60) %>%
bind_rows(sim_summary_70) %>%
bind_rows(sim_summary_80) %>%
bind_rows(sim_summary_100) %>%
bind_rows(sim_summary_150) %>%
bind_rows(sim_summary_365) %>%
mutate(immunity_duration = 1/omega) %>%
arrange(immunity_duration) %>%
mutate(model="single")
write_csv(waning_results_single, file = "/Users/matthewhoyle/Github_R_projects/Hunter_Gatherer_models/waning_results_single.csv")
waning_results_single
NA
ggplot(waning_results_single, aes(immunity_duration, sum_persist)) +
geom_line()+
geom_point()+
theme_bw()
###Set-up
# Define Paramenters
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Patch size
U <- length(patchPopSize) # Number of patches
initial_infected <- as.vector(rmultinom(1, 1, rep(0.5, U))) # Initial infected (initial infected patch randomly generated)
initial_infected_patch <- which(initial_infected > 0)
simName <- "SIRS metapopulation model" # Simulation name
tf <- 365*3 # Final time
# Agta Hunter-Gatherer contact rates
within_pop_contact = 1
between_pop_contact = 0.5/U # normalised by number of patches
#Create the named initial state vector for the U-patch system.
x0_meta <- unlist(lapply(
seq_len(U),
function(i){
c(patchPopSize[i] - initial_infected[i], initial_infected[i], 0, 0, patchPopSize[i])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(i) paste0(c("S","E","I", "R", "N"), i)))
# Define the state change matrix for a single patch
nu_meta <- matrix(c( -1, 0, 0, +1, +1, -1, 0, 0, 0, 0, # S
+1, -1, 0, 0, 0, 0, -1, 0, 0, 0, # E
0, +1, -1, 0, 0, 0, 0, -1, 0, -1, # I
0, 0, +1, -1, 0, 0, 0, 0, -1, 0, # R
0, 0, 0, 0, +1, -1 ,-1, -1, -1, -1), # N
nrow=5,byrow=TRUE)
# Define propensity functions
# Mass-action
a_meta <-
unlist(lapply(
seq_len(U),
function(patch) {
i <- patch
patches <- 1:U
#j <- if (patch == 1) U else patch - 1
other_patches <- patches[-i]
patch_beta <- c()
for(k in (1:(U-1))){
patch_beta[k] = paste0("+(beta_", other_patches[k],i, "*I", other_patches[k], "/N", other_patches[k], ")*S", i)
}
c(
paste0("(beta_", i, i, "*I", i,"/N", i, ")*S",i, paste0(patch_beta, collapse="")), # Infection
paste0("sigma*E", i), # Becomes infecious
paste0("gamma*I", i), # Recovery from infection
paste0("omega*R", i), # Loss of immunity
paste0("mu*N", i), # Births
paste0("mu*S", i), # Deaths (S)
paste0("mu*E", i), # Deaths (E)
paste0("mu*I", i), # Deaths (I)
paste0("mu*R", i), # Deaths (R)
paste0("alpha*I", i) # Deaths from infection
)
}
))
Define functions for calculating R0 from next-generation matrix
# Calculate R0 from NGM
R0ngm <- function(nextgen_matrix) {
eigenvalues = eigen(nextgen_matrix, only.values = T)
R0 = max(abs(eigenvalues$values))
return(R0)
}
beta.ngm <- function(beta_matrix) {
eigenvalues = eigen(beta_matrix, only.values = T)
beta_ngm = max(abs(eigenvalues$values))
return(beta_ngm)
}
#Collect parameters
parms_meta <- list(
sigma = 0.175, # E to I rate
gamma = 0.2, # I to R rate
omega = 1/180, # R to S rate
mu = demo_sum$Birth_rate_daily_Mean, # Birth/death rate per person per day
alpha = 1/1000)
# Define transmission terms and populate next-generation matrix
beta <- 0.6
nextgen_matrix <- matrix(nrow = U, ncol = U, data = 0)
beta_matrix <- matrix(nrow = U, ncol = U, data = 0)
for(i in 1:U){
for(j in 1:U){
parms_meta[[paste0("beta_",i,i)]] = within_pop_contact*beta
nextgen_matrix[i,i] = within_pop_contact*beta*(1/parms_meta$gamma)
parms_meta[[paste0("beta_",j,i)]] = between_pop_contact*beta
nextgen_matrix[j,i] = between_pop_contact*beta*(1/parms_meta$gamma)
nextgen_matrix[i,j] = between_pop_contact*beta*(1/parms_meta$gamma)
parms_meta[[paste0("beta_",j,j)]] = within_pop_contact*beta
nextgen_matrix[j,j] = within_pop_contact*beta*(1/parms_meta$gamma)
beta_matrix[i,i] = within_pop_contact*beta
beta_matrix[j,i] = between_pop_contact*beta
beta_matrix[i,j] = between_pop_contact*beta
beta_matrix[j,j] = within_pop_contact*beta
}
parms_meta[[paste0("N", i)]] = patchPopSize[i]
}
# Run simulations with the Direct method
set.seed(4)
out_meta <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Plot
plot_data_meta <- out_meta$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta <- ggplot(data = plot_data_meta, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 1, scales = "free_y")+
labs(x="Time (Days)",
y="Number of Individuals",
colour="State")+
theme_bw()
plot_meta
ggsave(filename = "meta_plot.pdf",
plot = plot_meta,
device = "pdf",
width = 7,
height = 8,
path = "/Users/matthewhoyle/Github_R_projects/Hunter_Gatherer_models")
## Table showing extinction/transmission info for each patch
extinct_data_meta <- out_meta$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N")),
persist = case_when(state=="I" & count > 0 ~ T,
state=="I" & count == 0 ~ F)) %>%
drop_na() %>%
select(patch, count, persist)
extinct_data_meta
beta_meta <- beta.ngm(beta_matrix)
paste0("Beta for whole system = ", beta_meta)
[1] "Beta for whole system = 0.857142857142857"
R0_meta <- R0ngm(nextgen_matrix)
paste0("R0 = ", R0_meta)
[1] "R0 = 4.28571428571429"
paste0("Actual number of infecteds at end of sim = ", sum(extinct_data_meta$count))
[1] "Actual number of infecteds at end of sim = 0"
# Total number of infecteds at the end of sim across all patches
sim_endpoint_meta <- as_tibble(out_meta$data) %>%
slice_max(t) %>%
distinct()
paste0("Did simulation run reach final endpoint?")
[1] "Did simulation run reach final endpoint?"
if (sim_endpoint_meta$t >= tf) {
print("Yes")
} else {
print("No")}
[1] "Yes"
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta <- list()
sim_list_meta <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta <- out_100_meta$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta[[i]] <- sim_data_meta
}
sim_output_meta <- bind_rows(sim_list_meta)
# Summary table of endpoint data
sim_output_meta <- sim_output_meta %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta
# Make Summary Table of output
sim_summary_meta <- sim_output_meta %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100,
mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/180)
sim_summary_meta
#Collect parameters
parms_meta_0 <- parms_meta
parms_meta_0$omega <- 0
# Run simulations with the Direct method
set.seed(4)
out_meta_0 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_0,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_0 <- out_meta_0$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_0 <- ggplot(data = plot_data_meta_0, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_0
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_0 <- list()
sim_list_meta_0 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_0 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_0,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_0 <- out_100_meta_0$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_0[[i]] <- sim_data_meta_0
}
sim_output_meta_0 <- bind_rows(sim_list_meta_0)
# Summary table of endpoint data
sim_output_meta_0 <- sim_output_meta_0 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_0
# Make Summary Table of output
sim_summary_meta_0 <- sim_output_meta_0 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100,
mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 0)
sim_summary_meta_0
#Collect parameters
parms_meta_1 <- parms_meta
parms_meta_1$omega <- 1
# Run simulations with the Direct method
set.seed(4)
out_meta_1 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_1,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_1 <- out_meta_1$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_1 <- ggplot(data = plot_data_meta_1, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_1
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_1 <- list()
sim_list_meta_1 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_1 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_1,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_1 <- out_100_meta_1$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_1[[i]] <- sim_data_meta_1
}
sim_output_meta_1 <- bind_rows(sim_list_meta_1)
# Summary table of endpoint data
sim_output_meta_1 <- sim_output_meta_1 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_1
# Make Summary Table of output
sim_summary_meta_1 <- sim_output_meta_1 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1)
sim_summary_meta_1
#Collect parameters
parms_meta_3 <- parms_meta
parms_meta_3$omega <- 1/3
# Run simulations with the Direct method
set.seed(4)
out_meta_3 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_3,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_3 <- out_meta_3$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_3 <- ggplot(data = plot_data_meta_3, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_3
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_3 <- list()
sim_list_meta_3 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_3 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_3,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_3 <- out_100_meta_3$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_3[[i]] <- sim_data_meta_3
}
sim_output_meta_3 <- bind_rows(sim_list_meta_3)
# Summary table of endpoint data
sim_output_meta_3 <- sim_output_meta_3 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_3
# Make Summary Table of output
sim_summary_meta_3 <- sim_output_meta_3 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/3)
sim_summary_meta_3
#Collect parameters
parms_meta_7 <- parms_meta
parms_meta_7$omega <- 1/7
# Run simulations with the Direct method
set.seed(4)
out_meta_7 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_7,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_7 <- out_meta_7$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_7 <- ggplot(data = plot_data_meta_7, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_7
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_7 <- list()
sim_list_meta_7 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_7 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_7,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_7 <- out_100_meta_7$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_7[[i]] <- sim_data_meta_7
}
sim_output_meta_7 <- bind_rows(sim_list_meta_7)
# Summary table of endpoint data
sim_output_meta_7 <- sim_output_meta_7 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_7
# Make Summary Table of output
sim_summary_meta_7 <- sim_output_meta_7 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/7)
sim_summary_meta_7
#Collect parameters
parms_meta_10 <- parms_meta
parms_meta_10$omega <- 1/10
# Run simulations with the Direct method
set.seed(4)
out_meta_10 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_10,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_10 <- out_meta_10$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_10 <- ggplot(data = plot_data_meta_10, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_10
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_10 <- list()
sim_list_meta_10 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_10 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_10,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_10 <- out_100_meta_10$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_10[[i]] <- sim_data_meta_10
}
sim_output_meta_10 <- bind_rows(sim_list_meta_10)
# Summary table of endpoint data
sim_output_meta_10 <- sim_output_meta_10 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
# Make Summary Table of output
sim_summary_meta_10 <- sim_output_meta_10 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/14)
sim_summary_meta_10
#Collect parameters
parms_meta_20 <- parms_meta
parms_meta_20$omega <- 1/20
# Run simulations with the Direct method
set.seed(4)
out_meta_20 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_20,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_20 <- out_meta_20$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_20 <- ggplot(data = plot_data_meta_20, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_20
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_20 <- list()
sim_list_meta_20 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_20 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_20,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_20 <- out_100_meta_20$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_20[[i]] <- sim_data_meta_20
}
sim_output_meta_20 <- bind_rows(sim_list_meta_20)
# Summary table of endpoint data
sim_output_meta_20 <- sim_output_meta_20 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_20
# Make Summary Table of output
sim_summary_meta_20 <- sim_output_meta_20 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/20)
sim_summary_meta_20
#Collect parameters
parms_meta_30 <- parms_meta
parms_meta_30$omega <- 1/30
# Run simulations with the Direct method
set.seed(4)
out_meta_30 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_30,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_30 <- out_meta_30$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_30 <- ggplot(data = plot_data_meta_30, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_30
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_30 <- list()
sim_list_meta_30 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_30 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_30,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_30 <- out_100_meta_30$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_30[[i]] <- sim_data_meta_30
}
sim_output_meta_30 <- bind_rows(sim_list_meta_30)
# Summary table of endpoint data
sim_output_meta_30 <- sim_output_meta_30 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_30
# Make Summary Table of output
sim_summary_meta_30 <- sim_output_meta_30 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100,
mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/30)
sim_summary_meta_30
#Collect parameters
parms_meta_40 <- parms_meta
parms_meta_40$omega <- 1/40
# Run simulations with the Direct method
set.seed(4)
out_meta_40 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_40,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_40 <- out_meta_40$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_40 <- ggplot(data = plot_data_meta_40, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_40
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_40 <- list()
sim_list_meta_40 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_40 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_40,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_40 <- out_100_meta_40$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_40[[i]] <- sim_data_meta_40
}
sim_output_meta_40 <- bind_rows(sim_list_meta_40)
# Summary table of endpoint data
sim_output_meta_40 <- sim_output_meta_40 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_40
# Make Summary Table of output
sim_summary_meta_40 <- sim_output_meta_40 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/40)
sim_summary_meta_40
#Collect parameters
parms_meta_50 <- parms_meta
parms_meta_50$omega <- 1/50
# Run simulations with the Direct method
set.seed(4)
out_meta_50 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_50,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_50 <- out_meta_50$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_50 <- ggplot(data = plot_data_meta_50, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_50
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_50 <- list()
sim_list_meta_50 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_50 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_50,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_50 <- out_100_meta_50$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_50[[i]] <- sim_data_meta_50
}
sim_output_meta_50 <- bind_rows(sim_list_meta_50)
# Summary table of endpoint data
sim_output_meta_50 <- sim_output_meta_50 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_50
# Make Summary Table of output
sim_summary_meta_50 <- sim_output_meta_50 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/50)
sim_summary_meta_50
#Collect parameters
parms_meta_60 <- parms_meta
parms_meta_60$omega <- 1/60
# Run simulations with the Direct method
set.seed(4)
out_meta_60 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_60,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_60 <- out_meta_60$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_60 <- ggplot(data = plot_data_meta_60, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_60
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_60 <- list()
sim_list_meta_60 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_60 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_60,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_60 <- out_100_meta_60$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_60[[i]] <- sim_data_meta_60
}
sim_output_meta_60 <- bind_rows(sim_list_meta_60)
# Summary table of endpoint data
sim_output_meta_60 <- sim_output_meta_60 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_60
# Make Summary Table of output
sim_summary_meta_60 <- sim_output_meta_60 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/60)
sim_summary_meta_60
#Collect parameters
parms_meta_70 <- parms_meta
parms_meta_70$omega <- 1/70
# Run simulations with the Direct method
set.seed(4)
out_meta_70 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_70,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_70 <- out_meta_70$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_70 <- ggplot(data = plot_data_meta_70, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_70
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_70 <- list()
sim_list_meta_70 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_70 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_70,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_70 <- out_100_meta_70$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_70[[i]] <- sim_data_meta_70
}
sim_output_meta_70 <- bind_rows(sim_list_meta_70)
# Summary table of endpoint data
sim_output_meta_70 <- sim_output_meta_70 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_70
# Make Summary Table of output
sim_summary_meta_70 <- sim_output_meta_70 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/70)
sim_summary_meta_70
#Collect parameters
parms_meta_80 <- parms_meta
parms_meta_80$omega <- 1/80
# Run simulations with the Direct method
set.seed(4)
out_meta_80 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_80,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_80 <- out_meta_80$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_80 <- ggplot(data = plot_data_meta_80, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_80
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_80 <- list()
sim_list_meta_80 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_80 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_80,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_80 <- out_100_meta_80$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_80[[i]] <- sim_data_meta_80
}
sim_output_meta_80 <- bind_rows(sim_list_meta_80)
# Summary table of endpoint data
sim_output_meta_80 <- sim_output_meta_80 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_80
# Make Summary Table of output
sim_summary_meta_80 <- sim_output_meta_80 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/80)
sim_summary_meta_80
#Collect parameters
parms_meta_90 <- parms_meta
parms_meta_90$omega <- 1/90
# Run simulations with the Direct method
set.seed(4)
out_meta_90 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_90,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_90 <- out_meta_90$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_90 <- ggplot(data = plot_data_meta_90, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_90
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_90 <- list()
sim_list_meta_90 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_90 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_90,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_90 <- out_100_meta_90$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_90[[i]] <- sim_data_meta_90
}
sim_output_meta_90 <- bind_rows(sim_list_meta_90)
# Summary table of endpoint data
sim_output_meta_90 <- sim_output_meta_90 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_90
# Make Summary Table of output
sim_summary_meta_90 <- sim_output_meta_90 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/90)
sim_summary_meta_90
#Collect parameters
parms_meta_100 <- parms_meta
parms_meta_100$omega <- 1/100
# Run simulations with the Direct method
set.seed(4)
out_meta_100 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_100,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_100 <- out_meta_100$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_100 <- ggplot(data = plot_data_meta_100, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_100
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_100 <- list()
sim_list_meta_100 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_100 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_100,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_100 <- out_100_meta_100$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_100[[i]] <- sim_data_meta_100
}
sim_output_meta_100 <- bind_rows(sim_list_meta_100)
# Summary table of endpoint data
sim_output_meta_100 <- sim_output_meta_100 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_100
# Make Summary Table of output
sim_summary_meta_100 <- sim_output_meta_100 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/100)
sim_summary_meta_100
#Collect parameters
parms_meta_110 <- parms_meta
parms_meta_110$omega <- 1/110
# Run simulations with the Direct method
set.seed(4)
out_meta_110 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_110,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_110 <- out_meta_110$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_110 <- ggplot(data = plot_data_meta_110, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_110
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_110 <- list()
sim_list_meta_110 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_110 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_110,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_110 <- out_100_meta_110$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_110[[i]] <- sim_data_meta_110
}
sim_output_meta_110 <- bind_rows(sim_list_meta_110)
# Summary table of endpoint data
sim_output_meta_110 <- sim_output_meta_110 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_110
# Make Summary Table of output
sim_summary_meta_110 <- sim_output_meta_110 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/110)
sim_summary_meta_110
#Collect parameters
parms_meta_120 <- parms_meta
parms_meta_120$omega <- 1/120
# Run simulations with the Direct method
set.seed(4)
out_meta_120 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_120,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_120 <- out_meta_120$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_120 <- ggplot(data = plot_data_meta_120, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_120
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_120 <- list()
sim_list_meta_120 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_120 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_120,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_120 <- out_100_meta_120$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_120[[i]] <- sim_data_meta_120
}
sim_output_meta_120 <- bind_rows(sim_list_meta_120)
# Summary table of endpoint data
sim_output_meta_120 <- sim_output_meta_120 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_120
# Make Summary Table of output
sim_summary_meta_120 <- sim_output_meta_120 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/120)
sim_summary_meta_120
#Collect parameters
parms_meta_130 <- parms_meta
parms_meta_130$omega <- 1/130
# Run simulations with the Direct method
set.seed(4)
out_meta_130 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_130,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_130 <- out_meta_130$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_130 <- ggplot(data = plot_data_meta_130, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_130
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_130 <- list()
sim_list_meta_130 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_130 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_130,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_130 <- out_100_meta_130$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_130[[i]] <- sim_data_meta_130
}
sim_output_meta_130 <- bind_rows(sim_list_meta_130)
# Summary table of endpoint data
sim_output_meta_130 <- sim_output_meta_130 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_130
# Make Summary Table of output
sim_summary_meta_130 <- sim_output_meta_130 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/130)
sim_summary_meta_130
#Collect parameters
parms_meta_150 <- parms_meta
parms_meta_150$omega <- 1/150
# Run simulations with the Direct method
set.seed(4)
out_meta_150 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_150,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_150 <- out_meta_150$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_150 <- ggplot(data = plot_data_meta_150, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_150
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_150 <- list()
sim_list_meta_150 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_150 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_150,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_150 <- out_100_meta_150$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_150[[i]] <- sim_data_meta_150
}
sim_output_meta_150 <- bind_rows(sim_list_meta_150)
# Summary table of endpoint data
sim_output_meta_150 <- sim_output_meta_150 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_150
# Make Summary Table of output
sim_summary_meta_150 <- sim_output_meta_150 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/150)
sim_summary_meta_150
#Collect parameters
parms_meta_220 <- parms_meta
parms_meta_220$omega <- 1/220
# Run simulations with the Direct method
set.seed(4)
out_meta_220 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_220,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_220 <- out_meta_220$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_220 <- ggplot(data = plot_data_meta_220, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_220
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_220 <- list()
sim_list_meta_220 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_220 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_220,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_220 <- out_100_meta_220$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_220[[i]] <- sim_data_meta_220
}
sim_output_meta_220 <- bind_rows(sim_list_meta_220)
# Summary table of endpoint data
sim_output_meta_220 <- sim_output_meta_220 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_220
# Make Summary Table of output
sim_summary_meta_220 <- sim_output_meta_220 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/220)
sim_summary_meta_220
#Collect parameters
parms_meta_270 <- parms_meta
parms_meta_270$omega <- 1/270
# Run simulations with the Direct method
set.seed(4)
out_meta_270 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_270,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_270 <- out_meta_270$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_270 <- ggplot(data = plot_data_meta_270, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_270
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_270 <- list()
sim_list_meta_270 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_270 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_270,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_270 <- out_100_meta_270$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_270[[i]] <- sim_data_meta_270
}
sim_output_meta_270 <- bind_rows(sim_list_meta_270)
# Summary table of endpoint data
sim_output_meta_270 <- sim_output_meta_270 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_270
# Make Summary Table of output
sim_summary_meta_270 <- sim_output_meta_270 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/270)
sim_summary_meta_270
#Collect parameters
parms_meta_365 <- parms_meta
parms_meta_365$omega <- 1/365
# Run simulations with the Direct method
set.seed(4)
out_meta_365 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_365,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
## Extra Plots
plot_data_meta_365 <- out_meta_365$data %>%
as_tibble() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(state = factor(state, levels = c("S", "E", "I", "R", "N"))) #%>%
#filter(state != "N")
plot_meta_365 <- ggplot(data = plot_data_meta_365, aes(x=t, y=count, colour=state))+
geom_line()+
facet_wrap(~factor(patch, levels = unique(patch)) ,ncol = 3, scales = "free_y")+
labs(x="Time",
y="Frequency")+
theme_bw()
plot_meta_365
## Run multiple simulations and saving output
num_sims <- 1000
sim_list_meta_365 <- list()
sim_list_meta_365 <- vector("list", length = num_sims)
for (i in 1:num_sims){
set.seed(i)
patchPopSize <- sample(camps.data$camp_total, 7, replace = TRUE) # Sample different patch sizes for each sim
x0_meta <- unlist(
lapply(
seq_len(U),
function(x){
c(patchPopSize[x] - initial_infected[x], initial_infected[x], 0, 0, patchPopSize[x])
}
))
names(x0_meta) <- unlist(lapply(seq_len(U), function(x) paste0(c("S","E","I", "R", "N"), x)))
out_100_meta_365 <- ssa(
x0 = x0_meta,
a = a_meta,
nu = nu_meta,
parms = parms_meta_365,
tf = tf,
method = ssa.d(),
simName = simName,
verbose = FALSE,
consoleInterval = 1
)
# Extract Final time point from output data
sim_data_meta_365 <- out_100_meta_365$data %>%
as_tibble() %>%
slice_max(t) %>%
distinct() %>%
pivot_longer(!t, names_to = "ID", values_to = "count") %>%
separate(ID,
into = c("state", "patch"),
sep = "(?<=[A-Za-z])(?=[0-9])") %>%
pivot_wider(names_from = state, values_from = count) %>%
mutate(persist = case_when(I > 0 ~ T,
I == 0 ~ F),
sim = i) %>%
select(sim, patch, I, N, persist)
sim_list_meta_365[[i]] <- sim_data_meta_365
}
sim_output_meta_365 <- bind_rows(sim_list_meta_365)
# Summary table of endpoint data
sim_output_meta_365 <- sim_output_meta_365 %>%
group_by(sim) %>%
summarise(total_I = sum(I),
total_N = sum(N)) %>%
mutate(percent_persist = total_I/(total_N)*100,
persist = case_when(total_I > 0 ~ T,
total_I == 0 ~ F))
sim_output_meta_365
# Make Summary Table of output
sim_summary_meta_365 <- sim_output_meta_365 %>%
summarise(mean_infecteds = mean(total_I),
sum_persist = (sum(persist, na.rm = T)/num_sims)*100, mean_percent_infected = mean(percent_persist)) %>%
mutate(omega = 1/365)
sim_summary_meta_365
Single
waning_results <- sim_summary_meta %>%
bind_rows(sim_summary_meta_3) %>%
bind_rows(sim_summary_meta_7) %>%
bind_rows(sim_summary_meta_10) %>%
bind_rows(sim_summary_meta_20) %>%
bind_rows(sim_summary_meta_30) %>%
bind_rows(sim_summary_meta_40) %>%
bind_rows(sim_summary_meta_50) %>%
bind_rows(sim_summary_meta_60) %>%
bind_rows(sim_summary_meta_70) %>%
bind_rows(sim_summary_meta_80) %>%
bind_rows(sim_summary_meta_90) %>%
bind_rows(sim_summary_meta_100) %>%
bind_rows(sim_summary_meta_110) %>%
bind_rows(sim_summary_meta_120) %>%
bind_rows(sim_summary_meta_130) %>%
bind_rows(sim_summary_meta_150) %>%
bind_rows(sim_summary_meta_220) %>%
bind_rows(sim_summary_meta_270) %>%
bind_rows(sim_summary_meta_365) %>%
mutate(immunity_duration = 1/omega) %>%
arrange(immunity_duration) %>%
mutate(model = "meta")
write_csv(waning_results, file = "/Users/matthewhoyle/Github_R_projects/Hunter_Gatherer_models/waning_results.csv")
waning_results
NA
ggplot(waning_results, aes(immunity_duration, sum_persist)) +
geom_line()+
geom_point()+
theme_bw()
combined_waning <- waning_results %>%
bind_rows(waning_results_single)
combined_waning
write_csv(combined_waning, file = "/Users/matthewhoyle/Github_R_projects/Hunter_Gatherer_models/combined_waning_results.csv")
combined_plot <- ggplot(combined_waning, aes(immunity_duration, sum_persist, colour = model))+
geom_line()+
geom_point()+
geom_segment(x = -Inf, y = 50, xend = 91.5, yend = 50, linetype = "dashed", colour = "grey") +
geom_segment(x = 5, y = 50, xend = 5, yend = -Inf, linetype = "dashed", colour = "grey") +
geom_segment(x = 91.5, y = 50, xend = 91.5, yend = -Inf, linetype = "dashed", colour = "grey") +
labs(x = "Duration of immunity",
y = "Probability of persistence after 3 years (%)",
colour = "Model Type")+
scale_color_discrete(type = wes_palette("Darjeeling1", type = "discrete")[1:2],
labels = c("Metapopulation", "Single Population"))+
theme_bw()
combined_plot
ggsave(filename = "combined_plot.pdf",
plot = combined_plot,
device = "pdf",
width = 7,
height = 5,
path = "/Users/matthewhoyle/Github_R_projects/Hunter_Gatherer_models")